Simulation of time evolution with the MERA 
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We describe an algorithm to simulate time evolution using the Multi-scale Entanglement Renor- 
malization Ansatz (MERA) and test it by studying a critical Ising chain with periodic boundary 
conditions and with up to L ~ 10® quantum spins. The cost of a simulation, which scales as L log(L), 
is reduced to log(L) when the system is invariant under translations. By simulating an evolution in 
imaginary time, we compute the ground state of the system. The errors in the ground state energy 
display no evident dependence on the system size. The algorithm can be extended to lattice systems 
in higher spatial dimensions. 

PACS numbers: 03.67.-a,70.00.00,02.70.-c 



The role of numerical simulations in many branches 
of physics is becoming more and more fundamental as 
the complexity of the systems of interest increases 
Recently, the injection of quantum information concepts 
has opened up the possibility of significant improvements 
in our ability to simulate strongly correlated quantum 
many-body systems. The entanglement present in the 
many-body wave function has been identified as a ma- 
jor limiting factor in numerical simulations. Accord- 
ingly, a big effort has been made within the quantum 
information community to devise new simulation strate- 
des that, building upon the matrix product state (MPS) 

and the density matrix renormalization group 
include a careful description of entanglement (see e.g. 
Refs. 4, 5, &, J., &, J., |lO|). In particular, the notion of 
entanglement renormalization — the systematic removal 
of short-ranged entanglement in the system — has been 
put forward as a means to obtain an efficient real-space 
renormalization group (RG) transformation for quantum 
systems on a lattice [7]. Relatedly, the Multi-scale En- 
tanglement Renormalization Ansatz (MERA) has been 
proposed as a variational many-body wave function to 
describe ground states [3]. It has been already demon- 
strated that the MERA offers a particularly accurate 
and compact representation of critical and non-critical 
ground states in ID lattices f^. However, the existence 
of an efficient algorithm to systematically compute the 
MERA for ground states is not yet demonstrated. 

In this paper we present an algorithm, referred to as 
the t-MERA algorithm, to simulate time-evolution with 
the MERA. For simplicity, we describe and test the ap- 
proach in a one-dimensional system, namely a critical 
Ising chain with periodic boundary conditions. However, 
the algorithm can be readily generalized to lattice sys- 
tems in higher dimensions. The cost of simulating L spins 
in an inhomogeneous system scales as Llog(L). We ex- 
ploit translational invariance to further reduce this cost 
to log(i), allowing us to accurately address systems of 
up to 2^° « 10^ spins with very modest computational 



resources. 

The t-MERA algorithm is inspired in the time-evolving 
block decimation (TEBD) algorithm for MPS d]. As in 
the latter, the tensors in the ansatz are updated so as to 
account for the action of a two-body gate acting between 
two neighboring lattice sites. However, while in an MPS 
the update involves only the tensor immediately close to 
those sites and is performed with a simple singular value 
decomposition 4], in the case of a MERA the update is 
given by a more sophisticated optimization, defined by a 
fidelity maximization, as described below. 

We consider a many body quantum system composed 
of i = 2^+^ sites, each of them described by a lo- 
cal Hamiltonian and some nearest neighbor interactions 
n^,^+l : C*' Hi]. The global Hamiltonian is then 
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where k ~ \,L with i + 1 = 1 for periodic boundary 
conditions. Consider the ensemble of wave functions that 
can be described exactly via a given MERA structure 
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unitary operator {xx^ = 
X'X ^ 1), ilJ,«J = -L m /3,, isometry (ftf = 1), 
J2 l-^P = 1 and a, /?, a', = 1, . . . , min(m, d^"'"*"^), where i 
counts the tensor network level and j the position on a 
given level 0, • The MERA structure is represented 
in Figure [T] a. As shown in Q the tensors x can be 
interpreted as disentanglers. The parameter m is the di- 
mension of the projected space and is related with the 
number of state kept TOg// in the reduced density ma- 
trix of half system in the Density Matrix Renormaliza- 
tion Group algorithm Indeed, a simple calculation 
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FIG. 1: (Color online) a) MERA tensor structure with periodic boundary conditions and 1 = 3. b) Schematic representation 
of the fidelity ([Sjl. Black (red) tensors represents the tensors inside the causal cone of \tp) Blue tensors (at sides) are 

outside the causal cone, thus they are contracted for free (c)). d-e) Schematic representation of the expressions (j4| and ([5|l. 



can show that the two quantities are related through 
mo// — m^^^°^^^~^\ Notice that the logarithmic scal- 
ing allows to describe th e q uantum correlations present 
in ID critical chains 13, 1^ - 



The simulation of time evolution is achieved by updat- 
ing the MERA when the operator U = e~*^* is applied 
to the state. When the time is real, i G i?, U is an 
unitary evolution, whereas when t is imaginary, t = —irj 
with rj Q R, U = e~^^ is an euclidean time operator that 
in the large t limit projects onto the ground state of H. 
We use the Suzuki- Trotter decomposition at a given or- 
der to obtain a sequence of two sites operators Uk,k+i to 
be applied to the initial ansatz ,15|. The problem is 
then reduced to the application of the operators Uk,k+i 
to the MERA and to absorb it, recovering the original 
structure with the minimum error. In other words, given 
a £ A4 and Uk.k+i we want to find \tp) E A4 such 
that 



max, T,^ 
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To perform such maximization efficiently, we implement 
a recursive procedure andan optimization of every ele- 
ment belonging to the causal cone The maximiza- 
tion is carried out for each tensor separately as no ex- 
act method is known (notice that in the TEBD algo- 
rithm this is performed exactly via a single SVD [J]). 
We first set lip) = \^) and we compute the trace fidelity 
^ = Tr{Uk,k+iP^'^) with p^'^ = as represented in 

Figure [T]b. As shown in [l6| the trace fidelity is equal to 
the contraction of the part of the tensor network inside 
the causal cone, i.e. the part of the network that can be 
influenced by a local operation which, at each level, is 
composed of maximum two unitaries x or three isome- 
trics r (see Fig. [1] c). Indeed, as shown in the figure, 
due to the properties of the tensor involved, the part of 
the network outside the causal cone is contracted for free. 
Thus, to compute !F only the reduced density matrix p^'^ 
of the two involved sites is needed and the application of 
any local operator Uk,k+i will result in a modification of 
the tensors inside its causal cone. Notice that the fidelity 



can be expressed as a function of the reduced density ma- 
trix at the upper level writing explicitly its dependence 
with the tensors at the last level, that is 

T = Tr{Uk,k+ix[k. i]x[k + 1, f]p%x[k. + 1, ^]^), 

(4) 

as shown in Figure [l]d. We first update a single tensor, 
e.g. x[k,l]'^ , contracting all the other tensors obtaining 



(5) 



as shown in Fig.[T]e [1^]. The maximum of the fidelity is 
given by x[^,^] = ^ where V is the unitary part of the 
polar decomposition of the matrix B = VP [1^] ■ We then 
write the analogous relation ([5]) related to the second ten- 
sor to be maximized and update x[k + 1,^]. The max- 
imization can be repeated until convergence is reached. 
We can now express explicitly the relation @ with its 
dependence from the isometrics of the last level Fs and 
the updated xs 



T = rr(J7fc,fc+iF[fc, e\T[k + 1, e\T[k -f 2, i] 

p};c^f[k,e\^f[k + i,e\^f[k + 2,e\^). 



(6) 



with Uk,k+i = x[k,i\^x[k + lJ]Wk,k+ix[k4]x[k+lJ]- 
To perform the fidelity maximization we repeat the pre- 
vious operations optimizing the Fs separately, defining 
at every optimization a new operator B and performing 
its polar decomposition. Notice that in this case B is 
a rectangular matrix and V is an isometry. Finally the 
same procedure is repeated for every level of the tensor 
structure, until the top is reached. Particular attention 
is needed to perform the update of the uppermost Fs and 
the vector A. In this case one can again write the prob- 
lem to be optimized in the same way as before. However 
both Fs can be updated together: once computed the 
operator B it can be Schmidt decomposed in two dif- 
ferent tensors each one defining the new F[j, l]s. The 
singular values obtained define the new vector A. In the 
case of euclidean evolution, a renormalization (enforcing 
Si I'^iP = 1) will take in account of the loss of norm of 




FIG. 2: (Color online) Convergence of the computed ground 
state energy E{L) for L = 2*, ^ = 8, 10 . . . 20 to the exact 
ground state energy Eex{L) with periodic boundary condi- 
tions Si. 



■0') due to the non unitarity of the euchdean operator 
U. We have build an approximation of the wave function 
Uk,k+i\4') which maximizes the fidelity J- and \ip) € M.. 
We then repeat the operations previously described for 
every operator Uk,k+i and for every Trotter step At until 
the desired convergence is reached. 

In conclusion, we sketch the algorithm scheme for the 
sake of clarity: 

1. Decompose the evolution operator in operators that 
act on nearest neighbor physical sites Uk,k+i with 
the Trotter decomposition at desired order. 

2. For every Uk,k+i and for every Trotter step, set 
IV'} ^ IV-") and perform the following actions: 

For every level i > 1 

(a) Find the optimal xs obtained by maximiza- 
tion of the fidelity ([4]). Replace the old xs 
with the new ones, xs, obtained as unitary 
part of the polar decomposition of the oper- 
ator B. If necessary repeat the maximization 
until convergence is reached. 

(b) Repeat step (a) to update the Fs maximizing 
the fidelity ([6]). Again repeat the process if 
needed. 

when the uppermost level is reached (i = 1) 

(c) Find the new isometrics and the new norm 
vector via a Schmidt decomposition of unitary 
part of the operator B. If euclidean evolution 
is simulated, renormalize the vector of the sin- 
gular values to obtain A. Set 1-0) IV")- 

3. When desired, perform one and two sites observ- 
ables measurement computing Tr{Ok,k' Pqc) as de- 
scribed in lldl. 
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FIG. 3: (color online) Energy E{V) as a function of the system 
size (L = 2*, £ = 8, . . . , 20) for the critical Ising model (blue 
circles), m = 4, dt = 0.1 and T; = 800. Red full line rep- 
resents the exact thermodinamical limit E'^ = — 4/7r. Inset: 
Difference of the computed energy with respect to thermodi- 
namical limit Ai5 = EiV) — E^ . The red line represents the 
exact scaling. 



The t-MERA algorithm, as described before, requires an 
impressive limited amount of resources (e.g. few hours 
on a laptop for L — 2^^ for the translational invariant 
case), polynomial both in memory and time as a func- 
tion of the system size and of the desired m: given L 
and m the memory resources scale as 0{m^LlogL) that 
is the memory needed to store the tensor structures [7|. 
The computational times are dominated by the tensor 
contractions needed to compute the operator B during 
the trace fidelity maximization, made of at most op- 
erations [iBl . As for every Trotter step one needs to com- 
pute 0{LlogL) different B operators (one for each link 
and level), the algorithms scales as 0{mPL log L). Even 
though the scaling with the projected size m is poly- 
nomial, it might still need a huge amount of computa- 
tional time for big m due to the high polynomial scaling 
power which might result in a limitation of this algorithm 
usefulness. However, as we show later on, already with 
TO = 4 we obtained very high precisions. More impor- 
tant and differently from previous proposed algorithms, 
the projected space size m needed to keep the error con- 
stant does not depend on the system size L, allowing to 
increase the size of the system up to thousands sites with 
only a linear-log cost in term of computational resources. 
Finally, if the Hamiltonian is traslational invariant, one 
can take advantage of this symmetry reducing the sim- 
ulation cost to 0(log L) allowing impressive system size 
to be studied as shown in the following. 

Results: We now apply the t-MERA algorithm to the 
study of the Ising chain ground state as a benchmark of 
the precision of the results that can be obtained. The 



Ising model is defined as 



W=-^/ia,^ + a,X+i, (7) 
<fc> 

where < k > describes periodic boundary conditions. 
The model is known to be critical for h — 1 and it can 
be solved exactly via the usual mapping to the fermionic 
operators f20j. In Fig.[2]we plot some typical convergence 
of the ground state energy as a function of the imaginary 
time T for the critical {h — 1) Ising model with peri- 
odic boundary conditions starting from the completely 
polarized state \ipo) = | | . . . f). We start with disentan- 
glers set to the identity and after some time Ti we switch 
on the optimization procedure for the xs as can be seen 
in Fig. [2] where a non smooth behavior is present. In 
Figure [3] we plot the finite size scaling energy resulting 
from our simulations at given final time Tf. the result- 
ing energy is E{oo) — 1.273229... with an error with 
respect to the exact solution of AE = 10~^. As clearly 
seen in the inset of Fig. [3] with this algorithm we can 
accurately reproduce the finite size scaling of properties 
such as the energy. Finally in FigH] A we show the error 
6E = E[L) — Eex{L) at given time Tf of the computed 
energy E{L) with respect to the exact energy for finite 
size Eex{L) as a function of the system size L: the error 
appears to be size independent. This is the most no- 
ticeable feature of the t-MERA algorithm and it reflects 
the underlaying MERA tensor structure 0]. In Fig[3]B 
we show the exponential dependence of the error with 
respect to the cut dimension m for L = 32: Increasing 
m the error decrease exponentially while the resources 
needed by the algorithm (memory and CPU-time) scale 
polynomially. We mention that similar scaling of the er- 
rors have been recorded for local and nearest neighbor 
observables (errors of order 10"'^ at Tf = 800) and for 
different critical models (data not shown). 

In conclusion, the idea introduced here works effi- 
ciently on every tensor network which has a finite size 
causal cone, that is, we can apply this algorithm to the 
proposed extension of 2D MERA Q structures with no 
fundamental changes. Moreover, the extensions of the t- 
MERA algorithm to include long range interactions and 
to study open system are also possible following [1, 12H . 
The exponential suppression of the error increasing the 
cut dimension m in comparison to the polynomial scal- 
ing of the resources also in critical systems are features 
that if confirmed in higher dimensionality, candidate this 
algorithm for the study of critical systems unaffordable 
with different methods. The limitation of this algorithm 
arises from the computational times needed, however, the 
code parallelization can be easily implemeted (2^ . 

After the completion of this work MERA tensor struc- 
ture have been used to describe ground state in 2D lat- 
tice [i^l and topological order [3]. SM and MR thank 
V.Giovannetti and R. Fazio for very useful comments and 
encouragement, I.Latorre for interesting discussions, M. 




FIG. 4: (color online) A. Error as a function of the system 
size L of the computed energy with respect to the exact one 
SE = E{L) - Eex{L), m = A, dt = 0.1 and Tf = 800 
B. Error SE as a function of m for L = 32, dt = 0.1 and 
Tf = 800. The dashed green line is an exponentail fit. 
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